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Abstract 


This  document  provides  a  detailed  mathematical  description  of  several  detection  algorithms  that 
have  been  implemented  in  the  Surface  Moving  Target  Indicator  (SMTI)  processor  gmtipro2  and 
tested  on  the  actual  Synthetic  Aperture  Radar  (SAR)  data  from  RADARSAT-2  acquired  in  the 
Moving  Object  Detection  Experiment  (MODEX)  mode.  The  common  feature  of  the  described  al¬ 
gorithms  is  their  Constant  False  Alarm  Rate  (CFAR)  property.  All  of  the  detectors  arc  generalized 
with  respect  to  the  number  of  available  channels.  Thus,  they  arc  defined  both  for  Moving  Object  De¬ 
tection  Experiment  1  (MODEX1)  and  Moving  Object  Detection  Experiment  2  (MODEX2)  modes, 
which  have  two  physical  and  four  virtual  channels,  respectively.  All  of  the  detectors  work  in  the 
image  domain,  but  can  be  combined  with  Doppler  Centroid  (DC)  offset  focusing  in  order  to  inte¬ 
grate  all  signal  energy  and  to  improve  the  probability  of  detection.  DC  estimation,  as  one  of  the 
most  important  auxiliary  algorithms,  is  also  presented  in  the  form  of  mathematical  equations.  The 
emphasis  of  the  report  is  entirely  on  the  mathematical  formalism  behind  the  computer  code  that 
implements  these  algorithms.  Open  literature  references  are  provided  regarding  the  derivation  and 
analysis  of  these  algorithms. 


Resume 


Le  document  decrit  de  faqon  detaillee  les  equations  qui  sous-tendent  plusieurs  algorithmes  de  de¬ 
tection  mis  en  ceuvre  dans  le  processeur  gmtipro2  indicateur  de  cibles  terrestres  mobiles  (ICTM)  et 
mis  a  l’essai  a  Taide  des  donnees  reelles  du  radar  a  synthese  d’ouverture  (SAR)  de  RADARSAT-2 
acquises  en  mode  MODEX  (Essais  de  detection  d’objets  en  mouvement).  Ces  algorithmes  parta- 
gent  une  caracteristique  :  leur  taux  de  fausse  alarme  constant  (TFAC).  Tous  les  detecteurs  peuvent 
traiter  sont  generalises  a  1’ ensemble  des  canaux  disponibles ;  ils  peuvent  done  traiter  les  donnees 
des  modes  MODEX  1  (deux  canaux  physiques)  et  MODEX2  (quatre  canaux  virtuels).  Tous  les  de¬ 
tecteurs  travaillent  a  partir  des  representations  en  images ;  cependant,  on  peut  les  combiner  a  la 
focalisation  du  decalage  du  centrolde  Doppler  (DC)  afin  d’integrer  toute  Tenergie  des  signaux  et 
ainsi  ameliorer  les  chances  de  detecter  toute  cible  mobile.  Comme  1’ estimation  du  DC  est  Tun  des 
algorithmes  auxiliaires  les  plus  importants,  il  est  aussi  presente  sous  forme  d’ equations.  Le  rapport 
est  entierement  consacre  au  formalisme  mathematique  qui  sous-tend  le  code  informatique  mettant 
en  ceuvre  ces  algorithmes.  Des  documents  de  source  ouverte  sont  cites  sur  la  derivation  et  1’ analyse 
de  ces  algorithmes. 
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Executive  summary 


A  Class  of  CFAR  Detectors  Implemented  in  the  SAR-GMTI 
Processor  gmtipro2  :  Mathematical  Formulation  of  the  Algorithms 

M.  Dragosevic,  W.  Burwash;  DRDC  Ottawa  CR  2013-089;  Defence  Research  and 
Development  Canada;  February  2015. 

Background:  A  modular  software  architecture  was  developed  at  Defence  Research  and  Develop¬ 
ment  Canada  (DRDC)  to  serve  as  a  robust  and  efficient  framework  for  a  variety  of  signal  processing 
algorithms.  The  intended  purpose  of  this  software,  known  as  gmtipro2,  is  both  operational  and 
research  oriented.  Besides  many  general-purpose  signal  processing  techniques,  this  architecture 
supports  a  number  of  detection  and  estimation  methods  that  were  developed  specifically  for  the 
RADARSAT-2  Moving  Object  Detection  Experiment  (MODEX)  data.  The  implementation  of  the 
detection  and  estimation  algorithms  started  from  the  most  basic  ones  and  progressed  tow  ards  more 
advanced  approaches.  For  research  purposes,  it  is  highly  desirable  to  provide  a  precise  mathematical 
description  of  the  implemented  detection  and  estimation  algorithms. 

Principal  results:  A  number  of  detection  algorithms  arc  presented  in  a  unified  and  consistent  way 
in  full  mathematical  detail.  The  key  equations  arc  provided  in  a  systematic  way  without  the  burden 
of  implementation  technicalities. 

Significance  of  results:  Detectors  considered  in  this  report  come  from  a  development  line  that  was 
not  previously  specified.  All  of  the  presented  algorithms  arc  the  result  of  incremental  upgrades  by 
the  authors  of  the  gmtipro2  softw  are.  As  such,  this  development  line  needed  to  be  well  documented. 
This  presentation  clarifies  the  main  underlying  concepts  including  parametric  and  non-parametric 
detection,  reduced  rank  (projection)  and  full  rank  processing,  selection  of  detection  thresholds, 
Doppler  offset  processing  and  Doppler  estimation.  Mathematical  expressions  presented  herein  com¬ 
plement  the  other  documentation  on  gmtipro2. 

Future  work:  A  comparison  of  different  algorithms  will  follow  as  a  separate  study.  It  will  compare 
the  algorithms  within  and  beyond  the  presented  class  using  simulated  and  real-world  moving  targets. 
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Sommaire 


A  Class  of  CFAR  Detectors  Implemented  in  the  SAR-GMTI 
Processor  gmtipro2  :  Mathematical  Formulation  of  the  Algorithms 

M.  Dragosevic,  W.  Burwash  ;  DRDC  Ottawa  CR  2013-089  ;  Recherche  et  developpement 
pour  la  defense  Canada ;  fevrier  201 5. 

Contexte  :  Recherche  et  developpement  pour  la  defense  Canada  (RDDC)  a  developpe  une  ar¬ 
chitecture  logicielle  modulaire  afin  de  creer  un  cadre  d’ utilisation  robuste  et  efficient  de  divers 
algorithmes  de  traitement  des  signaux.  Ce  logiciel,  nomrne  gmtipro2,  vise  a  servir  tant  a  des  fins 
operationnelles  qu’a  la  recherche.  En  plus  de  nombreuses  techniques  generates  de  traitement  des 
signaux,  cette  architecture  prend  en  charge  des  methodes  de  detection  et  d’estimation  propres  aux 
donnees  MODEX  de  RADARSAT-2.  Pour  mettre  en  ceuvre  les  algorithmes  de  detection  et  d’esti¬ 
mation,  on  a  commence  par  les  plus  simples,  puis  on  a  graduellement  mis  en  ceuvre  des  algorithmes 
de  plus  en  plus  complexes.  II  est  imperatif,  aux  fins  de  recherche,  de  decrire  fidelement  les  equations 
mathematiques  utilisees  dans  la  mise  en  ceuvre  des  algorithmes  de  detection  et  d’estimation. 

Principaux  resultats  :  Divers  algorithmes  de  detection  sont  presentes  de  faqon  uniforme  et  unifiee 
et  expliques  en  detail  sous  forme  mathematique.  Les  equations  les  plus  importantes  sont  systema- 
tiquement  enoncees  sans  se  soucier  des  details  de  mise  en  ceuvre. 

Portee  des  resultats  :  Les  algorithmes  de  detection  traites  dans  le  rapport  decoulent  d’efforts  de 
developpement  n’ayant  pas  ete  decrits  auparavant.  Cornme  tous  les  algorithmes  presentes  sont  le 
fruit  d' ameliorations  graduelles  par  les  createurs  du  logiciel  gmtipro2,  il  est  necessaire  de  bien 
decrire  les  resultats  de  ces  efforts.  Le  rapport  decrit  clairement  les  principaux  concepts  utilises,  no- 
tamment  la  detection  parametrique  et  non  parametrique,  le  traitement  de  rang  reduit  (projection)  et 
de  plein  rang,  le  choix  des  seuils  de  detection  ainsi  que  le  traitement  du  decalage  Doppler  et  l’esti- 
mation  Doppler.  Les  expressions  mathematiques  presentees  sont  un  complement  a  la  documentation 
de  gmtipro2. 

Recherches  futures  :  Une  etude  distincte  comparera  les  divers  algorithmes  presentes  dans  le  rap¬ 
port  et  d’autres  algorithmes  a  l’aide  de  cibles  mobiles  reelles  et  simulees. 
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1  Background,  Motivation  and  Objective 


Defence  Research  and  Development  Canada  (DRDC)  has  developed  a  powerful  and  versatile  Syn¬ 
thetic  Aperture  Radar  (S  AR)  Surface  Moving  Target  Indicator  (SMTI)  processor  known  as  gmtipro2. 
The  software  architecture  of  the  gmtipro2  processor  allows  the  user  to  build  a  variety  of  algorithms 
using  as  building  blocks  the  independently  developed  components,  modules  and  plug-ins.  An  algo¬ 
rithm  is  implemented  as  a  streamline  of  modules  with  their  dependent  pre-module  and  post-module 
plug-ins,  as  specified  by  the  meta-data  /  configuration  extensible  Markup  Language  (XML)  file. 
The  processor  can  be  used  for  innovative,  research  oriented,  as  well  as  standard  processing,  defined 
by  several  standard  XML  files. 

Efficient  usage  of  the  processor  depends  on  good  understanding  of  its  currently  available  compo¬ 
nents.  It  has  been  noted  that  many  researchers  prefer  a  rigorous  mathematical  formulation  of  the 
implemented  algorithms  to  any  other  description,  such  as  block  diagrams,  flow  charts,  pseudo-code 
or  verbal.  This  is  especially  true  for  certain  types  of  algorithms  such  as  those  described  in  this 
document. 

The  class  of  detectors  considered  in  this  report  is  not  new.  Their  basic  concept,  rooted  in  the  Likeli¬ 
hood  Ratio  Test  (LRT),  Maximum  Likelihood  (ML)  and  Generalized  Likelihood  Ratio  Test  (GLRT) 
approach,  is  well  known  and  well  published.  However,  there  arc  many  possible  variations,  specific 
details  and  enhancements,  which  actually  make  the  difference.  The  objective  of  this  document  is  to 
describe  the  implemented  algorithms  more  precisely  using  mathematical  formalisms  and  to  provide 
a  reliable  single-source  reference  for  these  algorithms. 

The  set  of  equations  presented  herein  describe  the  implemented  algorithms  very  closely  as  their 
mathematical  equivalents.  They  do  not  capture  the  organization,  parallel  execution,  synchronization, 
interfaces,  storage  requirements,  error  handling  or  any  other  technical  issue  related  to  the  gmtipro2 
design.  Moreover,  they  arc  not  necessarily  transcribed  into  the  C  language  corresponding  to  a  single 
module  or  plug-in  code.  Only  the  mathematical  concepts  arc  captured. 
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2  Notation  and  Key  System  Parameters 


This  section  summarizes  the  notation  that  is  used  throughout. 

2.1  Common  Symbols 

By  convention,  bold  lower  case  symbols  are  used  for  vectors  and  bold  upper  case  symbols  arc  used 
for  matrices. 

The  following  notation  is  used  for  common  operations  and  mathematical  units: 

(•)*  Complex  conjugate 

(•)T  Transpose  of  a  vector  or  matrix 

{■)n  Hermitian  transpose  (transpose  and  complex  conjugate) 

[•]„  n- th  component  of  a  vector 

(•)  Averaging 

{■)x  Averaging  over  argument  x 

j  Imaginary  unit 

I  Unit  (identity)  matrix 

The  following  symbols  arc  used  for  some  of  the  most  common  system  parameters: 

P  Number  of  Moving  Object  Detection  Experiment  (MODEX)  channels  (2  or  4) 

X  Radar  wavelength 

/p  Pulse  Repetition  Frequency  (PRF) 

D(-)  Baseline 

The  following  symbols  arc  used  for  some  of  the  most  important  data-set  parameters: 

/c  Doppler  Centroid  (DC) 

Vs  Satellite  speed  in  the  Earth-Centered,  Earth-Fixed  (ECEF)  referent  system 
Vr  Radial  speed  of  the  target 

Signal  and  image  coordinates  arc  generally  represented  by  the  following  symbols: 

£  Slow  time 
x  Fast  time 

v  Cross-range,  along-track,  azimuth  index 
y  Slant  range,  down  range  index 

Signal  and  image  processing  makes  use  of  the  following  standard  symbols: 
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*(■) 

Signal  sample 

s(-) 

Signal  vector 

"*(•) 

Mask  sample 

a2 

Variance  of  the  signal  or  a  signal  component 

C 

Covariance  matrix 

Px 

Projection  matrix 

b(-) 

Steering  vector 

More  details  regarding  notation  arc  provided  in  the  following  subsections. 

2.2  Data  Samples 

It  is  convenient  to  represent  the  samples  from  P  channels  by  a  vector  in  the  /'-dimensional  complex 
space  Cp,  for  example: 


s(0,x,y) 

s(x,y)  =  : 

_  s{P-l,x,y) 
[s  {x,y)]p  =  s(p,x,y) 


eC 


(1) 

(2) 


where  p  G  0, . . .  ,P  —  1  is  the  channel  index,  x  is  the  azimuth  coordinate  and  y  is  the  slant  range 
coordinate.  In  words,  the  p-th  component  of  the  vector  s(x.y)  is  the  sample  s(p,x,y )  from  channel 
p.  The  same  symbols,  .v ( • )  and  s(-),  will  be  used  in  different  algorithms  to  indicate  the  samples  to 
which  the  algorithm  is  applied  (original  samples),  which  may  be  at  any  state  of  processing.  Where 
the  processing  state  (or  a  processing  parameter)  of  the  original  samples  matters,  it  will  be  verbally 
stated,  as  well  as  by  additional  symbols. 


2.3  Data  Masking 

Many  of  the  implemented  algorithms  accept  data  masking,  which  is  used  to  indicate  invalid  or 
undesired  data  at  some  stage  of  processing.  Examples  include  invalid  image  portions  (wedges)  after 
SAR  focusing,  land  or  water  surfaces  within  algorithms  specialized  for  a  certain  type  of  clutter, 
pixels  outside  a  region  of  interest,  pixels  occupied  by  targets  when  background  statistics  are  being 
estimated  and  so  on.  The  symbol  m(x,y )  will  be  used  for  the  mask  of  any  nature  to  indicate: 

.  .  f  1  for  valid  samples 

m(x,y)  |  q  for  invalid  samples 

Data  masking  is  implemented  simply  as  multiplication  m(x,y)s(x,y)  to  discriminate  between  valid 
and  invalid  or  undesired  samples. 
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2.4  Averaging 


Data  averaging  will  symbolically  be  represented  by  (•)  and,  as  a  rule,  will  include  masking  (with 
default  masking  m (x.y)  =  1  for  for  all  range  and  azimuth  bins).  For  example 


(s (x,y)) 


I.tI: vm(x,y)s(x,y) 
ExL ym(x,y) 


(4) 


denotes  averaging  in  two  dimensions  over  all  valid  pixels  s(x,y).  A  subscript  is  used  to  specify 
averaging  in  a  certain  dimension.  For  example 


{s(p,x,y))p  = 


Lpm(x,y)s(p,x,y ) 

P 


(5) 


denotes  averaging  over  all  available  channels. 


2.5  Baseline 

Many  algorithms  use  the  Along-Track  Interferometry  (ATI)  baseline  as  a  parameter.  For  multi¬ 
aperture  modes  (Moving  Object  Detection  Experiment  2  (MODEX2)),  there  arc  different  ways  to 
define  the  baselines.  The  following  channel  enumeration  (indexing  order)  is  adopted  for  the  two- 
way  phase  center  separations  Dip)  in  the  antenna  frame  of  reference: 

D(0)  =0 

D(p)  >  D(q )  for  p  >  q 

For  the  most  used  MODEX  modes, 

D(p)  ~  PD,  (8) 

corresponding  to  the  model  discussed  in  [1],  except  that  in  this  report  the  distances  are  measured 
from  the  trailing  (aft-most)  phase  center. 

2.6  Steering  Vectors 

A  steering  vector  b(|3)  is  a  vector  function  of  the  scalar  variable  [3  that  takes  the  following  form: 

exp(;'47t^p) 

b(P)  =  exp(j'47i^P)  •  (9) 

_  exp  (y 471  P ) 


(6) 

(V) 
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3  Data  Model  in  the  Image  Domain 


All  algorithms  described  in  this  document  arc  applied  in  the  image  domain.  However,  the  image 
may  be  formed  using  an  offset  in  the  DC,  relative  to  the  DC  value  that  is  appropriate  for  the  clutter 
and  stationary  targets  [2,  3], 

This  section  briefly  presents  the  underlying  data  model,  which  is  necessary  for  a  better  understand¬ 
ing  of  the  implemented  algorithms,  especially  the  algorithms  for  clutter  suppression.  Some,  though 
not  all,  of  the  implemented  algorithms  presume  the  covariance  model  provided  in  this  section. 

3.1  Sample  Model 

Considering  that  all  focusing  algorithms  create  artifacts  such  as  azimuth  shift  of  moving  targets  and 
PRF  ambiguities,  the  underlying  model  of  the  image  samples  includes  several  components: 

s(p,x,y)  =  sm(p,x,y)  +  sc(p,x,y)  +  sa(p,x,y)  +  e(p,x,y)  (10) 

Basically,  at  a  given  pixel  position,  the  multi-channel  signal  sample  is  modeled  as  a  superposition  of 
the  moving  target  signature  sm(p,x,y),  clutter  sc(p,x,y)  and  additive  noise  e(p.x,y),  where  moving 
target  and  clutter  components  originate  from  different  geographical  locations.  Due  to  spectral  alias¬ 
ing,  ambiguities  of  moving  or  stationary  targets  may  appeal-  in  the  image  at  the  specific  distances 
from  the  principal  target  signatures.  In  the  model  (10),  the  ambiguities  from  strong  stationary  (or 
moving)  targets  (often  originating  from  a  man-made  construction)  are  represented  by  a  separate 
component  of  the  signal,  sa(p,x,y).  Other  ambiguities,  which  are  always  present  to  some  extent, 
may  be  absorbed  in  the  term  sc(p,x,y).  For  simplicity  (and  in  practice),  the  number  of  signal  com¬ 
ponents  in  (10)  can  be  limited  to  two  or  three  at  any  one  particular  pixel  although  in  theory  their 
number  may  be  large  [4],  The  two  components  that  are  always  assumed  present  are  the  system  noise 
and  clutter.  They  are  present  in  all  pixels.  In  some  pixels  there  may  be  more  significant  contributors. 
For  example,  a  third  term  may  be  a  single  moving  target  superimposed  on  the  clutter  and  noise.  At 
some  image  locations,  a  single  strong  ambiguity  may  be  superimposed  on  the  clutter  and  noise.  It 
is  possible  that  multiple  overlaid  moving  targets  get  imaged  at  the  same  pixel  location.  It  is  also 
possible  that  a  mix  of  targets  and  ambiguities  occupy  the  same  pixel.  However,  such  cases  are  not 
very  likely  and  will  not  be  considered  in  further  text. 

Under  idealized  conditions,  each  term  of  (10)  can  be  modeled  as  follows. 

After  focusing  and  coregistration,  the  moving  target  term  can  be  expressed  as: 

sm(x,y)  =5mbm  =  Smb  (11) 

with  unknown  parameters:  amplitude  Sm  and  ATI  phase.  The  latter  is  directly  proportional  to  the 
ratio  of  the  target’s  radial  speed  Vr  and  the  antenna  speed  Vs.  The  amplitude  may  be  considered 
unknown  but  deterministic  or  it  may  be  considered  random,  presumably  with  a  known  distribution. 
Unless  the  mover  is  a  point  target,  its  amplitude  ,Sm  may  vary  across  the  surface  of  that  same  tar¬ 
get.  For  a  small  distributed  target,  its  ATI  phase  is  assumed  to  be  spatially  constant.  For  a  large 


DRDC  Ottawa  CR  2013-089 


5 


distributed  target  such  as  a  ship,  its  ATI  phase  is  assumed  to  be  gradually  variable  over  the  image 
space  that  it  occupies. 

Similarly,  the  stationary  clutter  (the  background  target  signature)  may  be  modeled  in  the  same  way, 
with  the  ATI  phase  equal  to  zero: 


sc(. x,y)  =  Schc  =  Scb(0)  (12) 

where  the  amplitude  Sc  is  assumed  unknown,  either  deterministic  or  random.  The  probability  density 
function  (pdf)  associated  with  the  stochastic  clutter  model  is  typically  not  Gaussian  due  to  the 
image  texture.  Great  simplifications  in  the  detector  derivation  can  be  made  if  the  clutter  pdf  can  be 
considered  as  multivariate  Gaussian. 

Ambiguities  sa(p,x,y )  can  be  modeled  like  moving  or  background  targets  in  the  form: 

sa(x,y)=5aba  (13) 

where  the  specific  phases  in  ba  arc  similar  to  the  mover  ATI  phases,  but  they  depend  on  the  system 
parameters  and  the  applied  DC  offset  (positive  or  negative).  In  particular,  the  ambiguity  ATI  phase  of 
a  stationary  target  is  a  combination  of  two  terms  [5],  The  first  is  ±2nD(p)fp/Vs,  which  depends  on 
how  close  the  system  parameters  arc  to  the  DPCA  condition  (which  is  D  =  Vs/fp).  The  second  paid 
is  either  0  or  ±71,  introduced  by  transmitter  or  receiver  time  multiplexing  for  virtual  multichannel 
modes  (MODEX2  with  P  =  4).  For  a  mover  ambiguity,  its  ATI  phase  additionally  depends  on  the 
mover  radial  speed.  For  the  same  target,  Sa  increases  as  more  of  target  energy  appeal's  in  the  spectral 
alias  (greater  DC  mismatch). 

The  above  simplified  model  describes  the  various  components  as  each  belonging  to  a  subspace  of 
Cp,  spanned  by  b(|3)  or  ba. 

The  additive  noise  in  the  model  (10)  is  a  white  Gaussian  process  with  zero  mean  and  variance  G”. 
Noise  components  of  different  channels  are  statistically  independent  and  identically  distributed. 

3.2  Covariance  Matrix  Model 

If  a  stochastic  target  model  is  used,  the  target  covariance  matrix  defines  its  second  order  statistics. 
It  is  assumed  to  be  structured  as  follows: 


Cm  =  G^bmb"  (14) 

Detector  derivation  for  stochastic  targets  is  simplified  if  it  can  be  assumed  that  the  pdf  of  a  stochastic 
target  is  multivariate  Gaussian  with  zero  mean  and  the  structured  covariance  matrix  as  shown  in 
(14).  Then  the  pdf  depends  only  on  the  unknown  variance  c,n  and  the  unknown  normalized  radial 
speed  (as  a  parameter  in  bm). 

Similarly,  the  clutter  covariance  matrix,  used  in  the  stochastic  clutter  model,  may  be  structured  as 
follows: 

Cc  =  G^bcbf  (15) 
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Unlike  (14),  this  expression  has  only  one  unknown,  the  clutter  variance  G£.  This  simplified  struc¬ 
tured  model  is  exploited  in  some  algorithms  for  clutter  suppression.  In  other  algorithms,  no  particu¬ 
lar  structure  of  Cc  is  presumed.  In  fact,  the  model  (15)  is  not  necessarily  accurate  because  the  clutter 
itself  may  contain  a  multitude  of  moving  scattering  facets  (vegetation,  ocean  surface),  mixed  con¬ 
tributions  from  azimuth  and  range  ambiguities,  and  the  effects  of  uncalibrated  channel  imbalance. 
All  these  factors  work  to  increase  the  rank  of  Cc  from  1  up  to  possibly  full  rank. 

In  analogy  with  (14)  and  (15),  one  may  write: 

Ca  =  aababf  (16) 

when  the  ambiguity  is  believed  to  come  from  a  random  target  and  the  only  unknown  parameter,  Ga, 
relates  to  the  aliased  energy.  This  value  varies  for  a  variable  DC  offset. 

Since  signal  components  in  (10)  arc  statistically  independent,  the  covariance  matrix  C  can  be  writ¬ 
ten  as  a  sum  of  two  or  more  component  covariance  matrices.  Assuming  all  kinds  of  contributors, 
the  covariance  matrix  is  modeled  as: 

C  =  Cm  +  Cc  +  Ca  +  G“I  (17) 
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4  Sample  Covariance  and  Coherence  Matrices 


Sample  covariance  and  sample  coherence  matrices  arc  used  in  the  definition  of  many  implemented 
algorithms.  They  arc  estimated  from  the  data  by  massive  averaging  (in  both  dimensions  taking 
into  account  masking),  which  makes  the  algorithms  data-adaptive.  The  data  arc  required  to  be  in 
the  image  domain  for  all  of  the  algorithms  described  in  this  report.  Because  of  this,  the  sample 
covariance  and  coherence  matrices  arc  estimated  in  that  domain. 


4.1  Sample  Covariance  Matrix 

The  covariance  matrix  C  G  CPxP  is  estimated  by  averaging  over  all  currently  available  data: 


C  =  (s(x,y)s(x,y)H)  = 


Ex ,ym{x,y) 


This  is  the  sample  covariance  matrix  as  opposed  to  the  mathematical  expectation  of  s(x,y)s(x,y)n . 
In  the  definition  of  certain  algorithms,  the  inverse,  C1,  or  the  inverted  square  root,  C  l/2,  is 
needed.  The  latter  is  defined  (and  computed)  as: 


-,/2  =  IV1/2 . " 


where  A*  are  the  eigenvalues  and  u/,  arc  the  eigenvectors  of  C. 


4.2  Sample  Coherence  Matrix  and  Channel  Balance 

The  coherence  matrix  H  G  CPxP  is  estimated  by  averaging  over  all  currently  available  data  and  by 
normalizing  the  result: 


_ {s{p,x,y)s{q,x,y)*)Xy _ 

/(x(p,x,y)s(p,x,y)*)XtyV/(s(q,x,y)s(q,x,y)*)XJ 


As  already  mentioned,  masking  is  implied  in  the  averaging  operation. 


Coherence  is  computed  as  a  complex  scalar  based  on  three  averaged  inner  products.  The  same  three 
averaged  inner  products,  which  form  (20),  are  also  used  for  simple  (scalar)  channel  calibration.  The 
phase  of  the  numerator  shows  the  bulk  phase  imbalance  and  is  countered  by  the  opposite  calibration 
phase  between  channels  p  and  q.  The  ratio  of  the  two  factors  in  the  denominator  provides  the 
necessary  level  adjustment  between  these  channels. 
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5  Projection  Matrix  for  Clutter  Cancellation 


The  projection  matrix  P  is  implemented  as  the  projector  to  the  orthogonal  complement  of  the 
“clutter  subspace”  assuming  that  the  clutter  is  spanned  by  the  columns  of  a  matrix  Bc : 

Px  =  I-Bc(BfBc)^1Bf  (21) 

The  purpose  of  the  projection  matrix  P 1  is  to  cancel  clutter  or,  at  least,  the  main  components  of  the 
clutter. 

The  clutter  subspace  can  be  defined  in  several  ways  (depending  on  the  model),  using  a  deterministic 
(geometry-based)  method  or  an  adaptive  (balancing)  method.  Supported  options,  specified  by  user 
at  run-time,  are: 


Bc  = 


< 


[-b(ftt)-] 

b(0) 

b(0)bspa 

b(0)b*pa 


V  t>bal 


numerical 
geometry,  baseband 

geometry,  baseband  and  upper  ambiguity 
geometry,  baseband  and  lower  ambiguity 
balance 


(22) 


where  b(-)  is  a  steering  vector  as  defined  in  (9)  and  bspa  and  bbai  are  special  steering  vectors 
defined  below.  It  is  noted  here  that  bspa  and  bbai  are  closely  related  to  the  steering  vectors  ba  and 
bc,  respectively,  which  introduced  in  section  3.2  for  modeling  the  covariance  matrix.  The  current 
version  of  the  software  does  not  support  Bc  of  rank  larger  than  2. 


5.1  Clutter  Subspace  Defined  Numerically 

In  the  first  option  in  (22),  it  is  assumed  that  the  clutter  is  spanned  by  steering  vectors,  b ( (3* ) ,  where 
(3/(  for  0  <  k  <  P  are  given  numerically  from  an  external  source.  These  arc  the  ATI  values  that  will 
be  notched  out. 


5.2  Basic  DPCA  Defined  by  Geometry 

The  next  option  is  the  simplest  case,  corresponding  to  the  basic  DPCA  algorithm.  In  this  case,  it 
is  assumed  that  all  clutter  samples  are  scaled  versions  of  b(0).  The  justification  for  this  choice 
is  simple.  Adopting  the  models  (12)  and  (15),  clutter  samples  arc  assumed  to  be  identical  in  all 
coregistered  channels.  Then  the  projection  matrix  P  =  P0  is  given  by: 

p»=I'Ro)^(0)b(0,b(0)''  (23) 

This  type  of  projection  should  be  applied  to  the  data  samples  after  channel  balancing  and  spatial 
coregistration.  In  this  case,  Pxs(x,y)  removes  the  mean  (computed  over  channels)  (s(p,x,y))p  = 
b(0)ffs (x,y)/P  from  each  component  of  s(x,y).  In  other  words,  a  notch  is  placed  at  ATI  phase  0 
and  should  notch  out  the  baseband  clutter.  The  argument  of  the  steering  vector  is  implicitly  taken 
to  be  0  based  on  geometry  (all  channels  see  the  entire  surface  through  the  same  geometry). 
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5.3  Clutter  Components  from  Baseband  and  One  Ambiguity 
Defined  by  Geometry 

The  next  two  options  in  (22)  arc  modifications  of  the  simple  DPCA  case,  where  it  is  again  assumed 
that  the  samples  have  been  balanced  and  spatially  coregistered,  but,  possibly,  processed  with  a  DC 
offset,  causing  stronger  clutter  ambiguities.  Using  geometry  calculations,  it  is  easy  to  derive  the 
phase  mismatch  of  the  clutter  ambiguities,  as  a  result  of  improper  interpolation  [5],  In  fact, 

bspa  =  b*a  (24) 

with  reference  to  (13)  and  considering  only  ambiguities  from  stationary  targets.  The  elements  of  ba, 
i.e.  bspa  in  (22)  have  been  discussed  in  section  3.1.  In  some  more  detail,  they  are  given  by: 

[bspa]  p  =  exP  ( j2tt/p  i^spa  (p))  =  exp(j27t8spa(p))  (25) 

where  i^spa  (/?)  is  the  slow  time  shift  necessary  to  achieve  spatial  coregistration  between  channel  p 
and  channel  0  and,  in  the  equivalent  formulation,  8spa(p)  is  the  normalized  shift  within  the  data  grid 
(in  terms  of  sample  spacing  Ds  =  Vs/fp).  As  discussed  in  section  3.2,  these  normalized  shifts  arc  cal¬ 
culated  by  combining  two  parameters.  The  first  one  is  the  ratio  between  the  phase  center  separation 
D(p)  and  the  slow  time  spacing  Ds.  The  second  one  is  the  relative  sampling  delay  between  channels, 
in  terms  of  Pulse  Repetition  Interval  (PRI),  in  the  time  multiplex  that  is  used  for  MODEX2,  which 
is  ±0.5  or,  in  other  words,  ±50%  of  the  PRI  (the  sign  depends  on  the  channel  order  on  file).  The 
values  of  the  relative  shifts  8spa(p)  arc  the  same  values  as  used  for  spatial  coregistration  if  this  oper¬ 
ation  is  done  explicitly  (not  adaptively).  Besides  the  baseband,  one  of  the  ambiguities  is  selected  for 
cancellation.  Which  ambiguity  is  selected  for  cancellation  depends  on  the  sign  of  the  DC  offset.  If 
DC  offset  is  positive  (for  capturing  vehicles  with  negative  radial  speed),  the  upper  ambiguity  should 
be  cancelled.  For  negative  DC  offset,  the  lower  ambiguity  should  be  cancelled.  The  implementation 
allows  other  options,  such  as  notching  both  clutter  ambiguities  (bspa  and  b*pa),  but  not  the  baseband 
clutter.  However,  this  option  is  not  very  useful. 

5.4  Clutter  Notch  Positioned  by  Channel  Balancing 

The  final  option  in  (22)  is  based  on  the  measured  average  (bulk)  phase  imbalance  between  channels, 
c[)bai (p)  =  ^(s*{p1x,y)s(0,x,y)},  as  mentioned  in  section  4.2.  This  value  is  expected  to  be  0  after 
channel  balancing.  In  other  words, 

bbai  =  bc*  (26) 

where  bc,  as  estimated  from  the  data,  may  slightly  deviate  from  the  deterministic  vector  b(0).  Typ¬ 
ically,  the  data-based  bc  is  close  to  the  fixed  b(0),  but  there  may  be  some  exceptions.  In  particular, 
some  processing  chains  include  focusing  in  combination  with  coregistration  (interpolation)  with  a 
DC  offset  (suitable  for  movers)  after  original  channel  balancing.  The  DC  offset  then  causes  clutter 
ambiguities  in  the  background  image,  with  ambiguity  ATI  phases  ±27t8spa(p),  as  already  discussed. 
This  disturbs  the  channel  balance.  When  the  DPCA  condition  is  nearly  met  (8spa(p)  close  to  an  inte¬ 
ger),  it  is  possible  that  a  single  notch  placed  in  between  0  ATI  phase  and  ±2jtSspa(p)  may  provide  a 
good  compromise  between  suppressing  the  baseband  clutter  and  the  ambiguous  clutter.  This  option 
is  based  on  measuring  the  phase  imbalance  between  channels 

<h>ai(p)  =  As*(P,x,y,fc)s(0,x,y;fc))  (27) 
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where  it  is  symbolically  indicated  that  processing  (focusing  with  coregistration)  has  been  carried 
out  with  a  selected  DC  parameter  fc.  The  notch  is  placed  adaptively  to  restore  the  channel  balance: 

[bbai]p  =  exp(;'<|>baiO))  (28) 
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6  Test  Statistic  for  Quadrature  Detector  - 
Nonparametric  Detection 


This  section  describes  a  class  of  implemented  quadrature  detectors  that  arc  defined  without  any 
reference  to  the  mover  parameters.  However,  DC  offset  processing  can  be  used  even  with  this  type 
of  nonparamteric  detection  to  increase  the  chance  of  capturing  all  of  the  mover  energy  [3]  when 
forming  the  test  statistic. 


The  input  signal  s(x,y)  is  assumed  to  be  spatially  coregistered.  The  test  statistic  t(x.y)  is  evaluated 
for  every  pixel  according  to  the  following  formula: 


p-t 


t(x,y)=  2  £ 

P=Po 


Iffipfay)!2 

(|bfsp(x,y)|2) 


(29) 


where  weights  hp,  p  =  Pq,.  . .  ,P  can  be  defined  in  different  ways.  Averaging  in  the  denominator  of 
(29)  goes  over  all  valid  samples  and  serves  as  a  normalization  factor. 


A  variation  of  this  concept  permits  the  use  of  different  DC  offset  values  for  a  different  set  of  weights 
bp,  so  that  data  arc  refocused  (P  —  Pq  times)  prior  to  applying  each  of  the  weights  and  then  the  cor¬ 
responding  results  arc  summed  together.  This  option  is  explicitly  indicated  in  the  following  modifi¬ 
cation  of  (29): 


p-i 


t(x,y)  =  2  £ 
P=Pq 


|bfep(*>y;/c(bp))|2 
(|bp  Sp(x,y;/C(bp))|2) 


(30) 


where  the  relation  between  /c  and  hp  remains  to  be  defined  by  the  user  (as  the  designer  of  the  over¬ 
all  algorithm).  While  (29)  is  implemented  in  one  of  the  gmtipro2  modules,  the  modification  (30) 
can  be  implemented  (at  run  time)  through  an  iteration  of  several  modules. 


6.1  Projection  Based  Quadrature  Detector  -  DPCA 


In  the  framework  of  the  DPCA  approach,  orthogonal  vectors  b/;  (p  =  Pq,.  . .  ,P)  are  chosen  to 
span  the  “signal”  subspace  (the  orthogonal  complement  of  the  Po-dimensional  clutter  subspace) 
and  s p(x,y)  in  (29)  stands  for  the  image  pixels  focused  with  specific  processing  parameters  (such 
as  the  processed  Doppler  bandwidth  or  spectral  windows).  There  is  a  subtle  difference  between 
two  implemented  forms:  2a£T2s(x,y)//PQ  s(x,y)  and  (30).  In  the  modified  version  (30),  samples 
sp(x,y,fc(b(p))  arc  additionally  assumed  to  be  focused  using  a  specific  DC  offset,  which  is  ad¬ 
justed  to  the  weights  b;).  Partial  detection  results,  corresponding  to  each  hp  set,  arc  also  available. 


For  P=2,  the  only  option  is: 

b‘  =b(7t)  =  _  _i 

since  the  background  is  assumed  col  I  incar  with  bo  =  b(0). 


(31) 


For  P  >  2,  there  arc  various  possibilities.  An  example  of  a  suitable  set  of  vectors  b/;  is  defined  by: 


[b  p\q  = 


1 

7? 


cxp(  j2npq/ P)  for  q  =  0. . . 


,  P  —  1  and  p  =  1 , . . . ,  P  —  1 


(32) 
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This  is  the  set  of  Fourier  vectors  excluding  bo  =  b(0).  Partial  results  from  this  valiant  of  the  DPCA 
detector  are  very  close  to  the  Velocity  Synthetic  Aperture  Radar  (VSAR)  concept  proposed  in  [6]. 
In  the  special  case  of  P  =  4,  (32)  becomes: 


bj  b2  b3 


1 

2 


1  1 

j  -1 
-1  1 

~j  -1 


(33) 


and  its  partial  outputs  arc  tuned  to  three  unambiguous  radial  speeds.  Since  hp  and  bp  _p  are  approx¬ 
imately  matched  to  radial  speeds  of  the  opposite  signs  (especially  for  a  uniform  linear  array  as  in 
(8)),  it  makes  sense  to  use  opposite  DC  offsets  for  focusing  sp  and  s/>  .p.  In  other  words,  the  applied 
DC  offsets  and  the  applied  weights  b;)  within  the  same  algorithm  iteration  should  be  simultaneously 
matched  to  either  positive  or  negative  radial  speeds. 


Another  suitable  option  specifically  for  P  =  4  is: 


bi  b2  b3 


1 

2 


1  1  1 

1  -1  -1 

-1  -1  1 

-1  1  -1 


(34) 


These  vectors  arc  not  directional  (do  not  favor  positive  or  negative  radial  speeds)  and  no  meaningful 
selection  of  the  DC  offsets  can  be  recommended.  It  is  noted  that  bi  provides  the  largest  gain  for 
slow  targets  (of  either  sign)  within  this  set. 


6.2  Covariance  Based  Quadrature  Detector  -  Full  Rank 
Decorrelator 


For  this  method,  Pq  =  0  and 


b 


o 


bP_i  1  =  cr1/2 


(35) 


This  special  case  of  the  test  statistic  can  be  represented  in  a  different  form.  Considering  that  in  this 
case 


C  !/2s(x,y)  = 


bpS(x,y) 


and  because 


(C  1/2s(x,y)s(x,y)ffC  1/2)  =  I 

it  follows  that  diagonal  elements  of  (37)  satisfy 

(bpS(x,y)s(x,y)ffbp)  =  {\b”s(x,y)\2)  =  1 

and  then,  the  test  statistic  (29)  can  be  rewritten  as: 

t(x,y)  =  s(x,y)HC~1s(x,y) 


(36) 


(37) 


(38) 


(39) 
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This  form  of  the  test  statistic,  which  tends  to  decorrelate  channels,  is  mathematically  the  same 
as  that  of  a  ship  detector  used  on  polarimetric  imagery  [7]  (in  which  case  it  tends  to  decorrelate 
polarimetric  channels). 

The  quadrature  detector  based  on  (35)  (or  (39))  has  proven  very  useful  for  ocean  applications.  It 
may  also  be  useful  as  a  preliminary  detector  in  a  more  complex  algorithm.  For  example,  it  can  be 
used  to  quickly  detect  and  mask  strong  signatures  of  large  ships,  which  leads  to  a  more  accurate 
measurement  of  the  clutter  statistics  needed  for  further  application  of  more  advanced  algorithms 
[8]. 
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7  Test  Statistic  for  Combined  Detection  and 
Estimation 


The  input  signal  s(x,y)  is  assumed  to  be  spatially  coregistered.  The  test  statistic  t(x.y)  is  evaluated 
for  every  pixel  according  to  the  following  set  of  formulae: 


J{x,y,vv-,fc)  = 


X^eiBy  I b(Vr/Vs)  Q(P  —  Pq\  fc)s(x  +  dx,y  +  dy\  fc)\" 
b{Vr/Vs)HQ(P  -  Po:fc)b(Vr/Vs) 


Vr(fc)  =  arg  max  J(x,y,VT;fc) 

VrGBy 

t(x,y)  =  max  (/(Vr(/C)  G  Dyr(/c))/(x,v,  t>r;/c)) 

/c6%c 


(40) 

(41) 

(42) 


In  accordance  and  in  addition  to  the  general  notation  introduced  in  Section  2,  the  following  symbols 
are  used: 


Wfe 

s(-;/c) 

vr 

Dv 

Bvr(/c) 

dx 

dy 

By 

/(•) 

b(-) 

Q(*;/c) 


The  interval  of  test  values  for  /c 
Image  pixel  focused  with  a  certain  fc  value 
Test  radial  speed  of  a  moving  target 
The  interval  of  unambiguous  values  for  Vr 

The  interval  of  test  values  for  Vr  to  be  acceptable  depending  on  the  fc  parameter 

Offset  in  the  neighborhood  of  the  current  azimuth 

Support  window  in  azimuth 

Offset  in  the  neighborhood  of  the  current  range 

Support  window  in  range 

Indicator  function  taking  value  1  for  true  arguments  and  0  for  false  arguments 
Test  steering  vector 

Kernel  matrix  of  rank  k,  may  vary  with  the  applied  fc 


The  kernel  Q (P  —  Po;fc),  used  to  suppress  the  clutter  in  (40),  is  redefined  for  each  of  the  specific 
methods  of  the  same  class.  It  may  have  full  rank  P  or  reduced  rank  P  —  Pq  where  Pa  is  the  assumed 
rank  of  the  clutter  subspace.  It  may  be  estimated  or  constructed  for  each  test  value  of  fc.  It  is 
emphasized  here  that  the  same  kernel,  which  is  used  for  clutter  suppression,  must  also  be  used  in 
normalization  (the  denominator  of  (40)). 


In  (40),  summing  the  pixel  contributions  within  a  rectangular  window,  centered  around  (jc,y),  and 
defined  by  Dx  and  Dv  is  optional  (set  through  design  parameters  at  run-time).  The  default  window  is 
one  sample  wide,  i.e.  Dx  =  {0}  and  Dv  =  {0}.  This  choice  (and  default  setting)  forces  the  algorithm 
to  use  a  single  pixel  value  s (x,y;fc)  when  evaluating  the  criterion  J(x,y,Vr;fc),  which  is  then  used 
to  estimate  the  radial  speed  in  (41)  and  to  compute  the  test  statistic  t(x,y )  in  (42). 


Focusing  may  be  done  with  a  variable  DC  offset,  while  the  other  focusing  parameters  (Doppler 
rate,  processed  bandwidth  and  spectral  window  shape)  arc  kept  constant  according  to  the  general 
algorithm  formulation  (40),  (41),  (42).  Normally,  these  fixed  parameters  take  the  appropriate  or 
default  values,  but  they  too  may  be  modified  by  the  user  as  design  parameter  at  run  time.  Multiple 
DC  offset  values  can  be  used  iteratively  in  order  to  evaluate  t(x.y)  defined  in  (42).  The  set  of  suitable 
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test  DC  offsets  is  provided  to  the  algorithm  externally  (as  a  design  parameter)  either  directly  as 
frequency  offsets  from  the  dominant  (bulk  clutter)  DC  or  as  a  set  of  preliminary  test  Vr  values. 
These  preliminary  test  values  arc  supposed  to  be  given  in  coarse  steps.  For  each  DC  offset  (or 
preliminary  test  VT),  a  large,  denser,  set  of  Vr  values  is  tested  on  a  finer  scale  (in  finer  steps)  using 
(41).  The  fine-scale  estimates  from  (41)  are  then  checked  against  the  expected  domain  D vr(fc)  to 
see  if  they  arc  actually  matching  the  coarse-level  setting  used  in  focusing.  If  this  is  the  case,  the 
estimate  produced  by  (41)  is  accepted,  otherwise  it  is  rejected.  In  the  latter  case,  it  is  very  likely 
that  the  non-matching  estimate  originates  from  a  mover  ambiguity.  By  this  method,  the  detection 
candidates  arc  censored.  Partial  maximization  results  (uncensored)  from  (41)  for  each  fixed  fc  arc 
also  available  as  algorithm  output. 

The  test  statistic  t(x,y )  of  (42)  is  used  for  moving  target  detection  by  comparison  to  a  suitable 
threshold.  The  value  Vr  that  achieves  the  maximum  globally  over  all  fc  is  an  estimate  of  the  target’s 
radial  speed. 

7.1  Reduced  Rank  Kernel  -  DPCA  with  Adaptive  Steering 

One  option  for  choosing  the  kernel  matrix  is  to  notch  out  the  clutter  using  any  one  form  of  the 
projection  matrix  (based  on  numerical  input,  geometry  calculations  or  balance  measurements)  as 
discussed  in  Section  5  and  then  to  normalize  it  by  an  estimate  of  the  residual  power  d2,  which  may 
change  with  /c. 

Q(P-Po)  =  2d~2Px  (43) 

The  reduced  rank  of  the  kernel,  P  —  Pq  <  P,  is  as  discussed  in  section  6.1.  Additionally,  the  projec¬ 
tion  matrix  may  change  with  fc  if  it  is  derived  from  the  balance  results,  (jbai  (p)-  If  the  projection 
matrix  is  based  only  on  the  geometry  considerations,  then  it  may  depend  on  the  sign  of  the  DC  off¬ 
set  (the  difference  between  fc  used  in  focusing  and  the  DC  value  corresponding  to  the  bulk  ground 
backscatter),  as  noted  in  Section  5. 

This  type  of  detector  is  described  in  detail  in  [9],  In  broad  terms,  it  is  closely  related  to  the  Notch 
Periodogram  (NP)  concept  [10,  11]. 


7.2  Full  Rank  Kernel  -  Adaptive  Matched  Steering 


The  option  given  by 


Q(p)  =  cr1 


(44) 


is  to  decorrelate  or  whiten  the  channels  using  the  sample  covariance  matrix  as  discussed  in  Section 
4.1.  With  this  option,  C  defined  in  Section  4.1  actually  changes  with  fc  due  to  a  change  in  the  ratio 
of  the  ambiguity  to  baseband  power. 


This  type  of  test  statistic  is  discussed  in  numerous  papers  (e.g.  [12]). 
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8  Detection  Thresholds 


Detection  is  declared  when  the  test  statistic  t(x,y)  surpasses  the  detection  threshold  Tj.  With  the 
cuiTent  implementation,  the  threshold  can  be  set  directly  or  indirectly. 

All  test  statistics  presented  in  sections  6  and  7  arc  normalized  to  make  the  choice  of  the  threshold 
easier.  If  r|  is  set  externally,  directly  as  a  number,  it  may  loosely  be  interpreted  as  the  relative  level 
above  the  power  of  the  residuals  after  clutter  suppression. 

A  higher  threshold,  used  on  the  same  test  statistic,  decreases  the  Probability  of  False  Alarms  (PFA), 
but  it  also  decreases  the  Probability  of  Detection  (PD).  The  relationship  between  r|  and  PFA  (or 
PD)  for  a  given  form  of  test  statistic  depends  on  the  pdf  of  the  measurements  and,  in  general, 
is  very  difficult  to  derive  analytically.  By  setting  the  allowed  PFA,  the  choice  of  r|  is  indirectly 
made  through  certain  implemented  expressions  that  arc  not  necessarily  accurate.  The  implemented 
expressions  arc  set  for  special  cases  and  under  special  conditions.  When  applied  to  other  cases,  they 
produce  a  value  for  r|  which  does  not  necessarily  guarantee  the  preselected  PFA. 

The  details  regarding  the  threshold  depend  on  the  type  of  the  test  statistic  and  differ  between  the 
non-parametric  and  parametric  detectors,  projection-based  and  full-rank  clutter  suppression,  single¬ 
pixel  and  multi-pixel  test  statistics.  The  common  assumption  for  all  is  the  assumption  of  Gaussian 
pdf  of  the  single  pixel  complex  residuals. 


8.1  Detection  Threshold  for  Nonparametric  CFAR  Detectors 

The  detection  threshold  r|  is  determined  to  satisfy  the  equation: 


1  -  Pfa  =  /  po(z;v)dz 
Jo 


(45) 


where  Pfa  is  the  chosen  probability  of  false  alarms,  v  =  2 n  is  the  known  or  estimated  value  of  the 
Degrees  of  Freedom  (DOF)  and  po(z)  is  the  assumed  pdf  of  the  statistic  z  under  the  hypothesis  that 
no  moving  target  is  present.  In  the  current  implementation,  it  is  assumed  that  this  pdf  is  %2,  which 
has  the  following  form 


po(z;v) 


zn~  Vz/2 

2  'T(n) 


(46) 


where  r(-)  is  the  gamma  function  and  the  DOF  depends  on  the  method  used  in  computing  the  statis¬ 
tic  z.  For  example,  if  the  projection-based  quadrature  detector  of  section  6.1  is  used,  then  z  =  t(x,y ) 
has  2(P  —  Pq)  DOF  (where,  typically  Po  =  1,  or  maybe  P(l  =  2).  For  this  detector,  the  test  statis¬ 
tic  is  designed  to  cancel  clutter  leaving  only  the  Gaussian  thermal  noise.  Under  such  assumptions 
of  Gaussian  samples  (in  the  absence  of  targets),  the  test  statistic  has  indeed  the  %2  pdf.  If  the  full 
rank  non-parametric  detector  of  section  6.2  is  used,  the  same  approach  is  valid  (with  p}  =  0)  if  the 
additional  assumption  is  made,  namely  that  the  residual  clutter  is  Gaussian. 


In  these  cases,  the  DOF  can  be  set  directly  or  estimated  from  the  data,  as  stated  in  section  8.3. 
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8.2  Detection  Threshold  for  Parametric  CFAR  Detectors  with 
Adaptive  Steering 

In  the  case  of  the  detector  described  in  Section  7,  taking  z  =  J{x,y,  Vr ;/c)  one  would  expect  the  DOF 
to  be  2  (for  no  spatial  averaging),  or  to  be  related  to  the  estimated  Equivalent  Number  of  Looks 
(ENL)  according  to  Section  8.3  (after  spatial  averaging).  The  %2  property  would  be  justified  either 
assuming  full  clutter  suppression  via  projection  or  assuming  that  the  residual  clutter  is  Gaussian. 
Either  way,  the  relationship  between  r|  and  Pfa,  given  by  (45)  and  (46)  is  good  only  for  the  parti al 
criterion  (40).  The  relationship  is  good  for  any,  arbitrary  value  of  Vr  in  (40),  but  it  is  not  valid  for 
the  final  test  statistic  (42). 

Nonetheless,  if  the  threshold  is  given  indirectly,  it  is  computed  as  in  section  8. 1  and  it  provides  the 
prescribed  Pfd  per  test  case,  but  the  overall  PFA,  /Ji  \  is  actually  higher. 


8.3  Equivalent  Number  of  Looks 


The  implemented  method  for  estimating  the  ENL,  n,  and  DOF,  v  =  2 n,  is: 


{S2(p,x,y)) 
S{p,x,y)  =  \s{p,x,y)\2 


where  averaging  and  masking  arc  applied  as  explained  in  Section  2. 


(47) 

(48) 
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9  Doppler  Estimation  and  Analysis 


The  DC  offset,  which  is  used  in  some  detectors,  is  applied  relative  to  the  DC  of  the  background 
clutter.  This  section  described  the  DC  estimator  that  is  implemented  to  estimate  the  bulk  clutter  DC. 

The  DC  estimator  uses  the  phase  increment  method  to  compute  the  localized  DC,  estimated  inde¬ 
pendently  at  different  fast  time  intervals.  Then,  an  over-all  DC  model  is  fitted  to  the  local  estimates. 


9.1  Local  DC  Estimation 

The  average  phase  increment  is  based  on  the  autocovariance: 

C(x)  =  ^m(C,x)m(C  +  AC,T)/(C,'tMC  +  A^T),  (49) 

? 

which  is  also  called  the  Average  Cross  Correlation  Coefficient  (ACCC)  [13].  Here,  All,  is  the  PRI. 
The  angle  (or  phase)  of  the  ACCC  is  given  by 

4*accc(x)  =  Z  (c(x))  .  (50) 


The  “raw”  baseband  DC  estimate  is  then  computed  by 

/c(x)  =  ^_4,accc('t)  (51) 

271 

9.2  Smooth  DC  Estimates 


A  suitably  biased  version  of  this  estimate  is  expressed  as: 

fi 


/b(x)  =  ^Z  (ejis> acccWg  7'ej 

(cOS^acccCr))^ 


cos(0) = 
sin(0)  = 


\J  (cos(<|)accc(x)))2  +  <sin(4>accc(x))) 
(sin(<|)accc  (x))^ 

\J  (cos(<|)accc(x))}^  +  <sin(4>accc(x))) 


(52) 

(53) 

(54) 


where  averaging  goes  over  all  processed  slow  time  £  and  implies  masking,  as  described  in  Section  2. 
The  introduced  bias  is  towards  0  and  away  from  ±/p/2,  which  serves  to  avoid  ambiguities  in  further 
processing  [14].  Linear  fitting  with  respect  to  fast  time  x  and  removal  of  outliers  arc  performed  on 
the  biased  sequence  /b(x)  defined  by  (52),  not  on  the  original  unbiased  sequence  /c(x)  of  (51).  The 
domain  of  valid  samples  Dx  is  defined  as: 


=  |x|(/b(x)  -/lin(x))2  <  a((/b(x)  -/lin(x))2)x}  (55) 
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where  fnn  (x)  is  the  linear  fit  to  the  biased  estimates  /b  (x)  as  defined  by  (52)  and  a  is  a  given  constant. 
The  sequence  of  outlier-free  DC  estimates,  /c(x),  is  then  obtained  as: 


/c(x)=/b(x)  +  ^0,  xeDt  (56) 

271 

and  the  final  fit  is  performed  on  /c(x)  over  the  domain  Dx  to  compute  the  coefficients  of  the  Doppler 
polynomial. 

For  the  computation  of  the  DC  map,  as  paid  of  DC  analysis,  averaging  is  applied  in  both  directions 
within  small  moving  windows: 

/c(x,C)  =  Z((^,x));(^  +  AC,x))Xf)Cc,  (57) 

where  /(x,  Q  is  the  Doppler  centroid  estimate  calculated  from  the  xr  x  i^c  data  chip  centered  at 

(x,Q- 
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10  Conclusion 


Detectors  described  in  this  report  represent  one  class  of  detectors  that  are  implemented  in  gmtipro2. 
They  have  many  commonalities.  They  all  include  a  mechanism  of  suppressing  the  clutter  based  on 
known  global  properties  of  the  clutter,  which  arc  given  a  priori  or  previously  estimated  on  the  entire 
data  set. 

Within  this  class,  the  described  detectors  can  be  divided  into  two  groups: 

-  DPCA  detectors  arc  based  on  matrix  projection;  they  make  use  of  the  assumption  that  the 
clutter  lies  in  the  clutter  subspace  of  a  low  rank  (1  or  2),  which  is,  as  a  rule,  a  priori  known. 

-  Full  rank  detectors  make  use  of  the  noisy  clutter  covariance  matrix,  which  needs  to  be  esti¬ 
mated  from  the  data  using  global  averaging. 

Another  classification  is  also  applicable  to  the  described  detectors: 

-  Nonparametric  detectors  arc  based  on  a  test  statistic  in  the  quadratic  form  without  any  refer¬ 
ence  to  any  signal  parameters. 

-  Parametric  detectors  imply  certain  target  parameters;  since  true  target  parameters  arc  a  priori 
unknown,  these  detectors  typically  perform  multiple  tests  and  maximization  over  the  param¬ 
eter  space. 

Yet  another  view  can  be  added  to  this  type  of  detectors: 

-  Single-sample  detectors  work  on  a  single  image  sample  to  generate  the  test  statistic,  which  is 
then  compared  to  a  threshold. 

-  Multiple-sample  detectors  can  use  several  neighboring  samples  to  compute  the  test  statistic; 
they  work  on  sliding  windows  and  their  thresholds  must  be  adjusted  accordingly. 

All  of  the  detectors,  by  their  design,  can  benefit  from  iterative  processing  with  different  DC  offsets. 

The  described  class  of  algorithms  offers  a  great  degree  of  flexibility  at  run  time.  To  exploit  the  full 
potential  of  the  implemented  gmtipro2  modules  and  plug-ins,  the  user  must  become  familial-  with 
their  functional  formulation,  but  not  necessarily  with  their  implementation  aspects. 

Evaluation  of  the  algorithm  performance  is  beyond  the  scope  of  this  document.  Two  examples  are 
shown  in  the  Annex  for  illustration  purposes  only.  The  first  example  illustrates  the  algorithm  effec¬ 
tiveness  for  land  applications,  and  the  second  demonstrates  the  algorithm  application  in  a  maritime 
scene.  Both  examples  show  the  results  of  combined  detection  and  estimation  of  moving  vehicles 
as  delivered  by  the  gmtipro2  processor.  The  estimated  target  locations  are  shown  overlaid  on  the 
background  image,  color  coded  according  to  their  estimated  radial  speeds. 
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Annex  A:  Algorithm  Output  Examples 


Two  processing  examples  are  shown  to  illustrate  the  output  of  the  presented  combined  detection 
and  estimation  algorithm.  Both  outputs  were  generated  using  the  gmtipro2  processor  on  selected 
RADARS  AT-2  data. 

Fig.  A.  1  is  a  composite  image  of  the  Ottawa  region.  The  background  is  an  optical  image  from  a 
commercial  source,  which  also  includes  the  roadmap.  The  overlaid  detections/estimations  of  mov¬ 
ing  vehicles  come  from  gmtipro2.  The  detection/estimation  results  arc  accumulated  from  multiple 
RADARS  AT-2  ascending  and  descending  passes  in  MODEX1  and  MODEX2  modes.  Detected  ve¬ 
hicles  are  presented  by  color-coded  points  where  color  represents  their  estimated  speed.  Traffic 
patterns  can  be  observed  in  a  qualitative  way.  This  is  a  reprint  of  the  poster  prepared  for  the  EUSAR 
2012  best  image  contest. 

Fig.  A.2  is  a  SAR  image  of  the  Gibraltar  region  with  overlaid  ship  detections/estimations.  Both 
the  SAR  image  and  the  ship  results  were  produced  by  gmtipro2  from  an  evening  RADARSAT-2 
pass  in  MODEX1  mode.  Most  ships  arc  visible  in  the  SAR  image.  They  are  shifted  in  azimuth 
from  their  estimated  positions.  Detected  vessels  arc  presented  by  color-coded  points  where  color 
represents  their  estimated  speed.  The  observed  azimuth  shifts  correspond  to  these  estimated  speeds 
in  sign  and  magnitude.  Littoral  and  ship  ambiguities  arc  also  visible  in  the  SAR  image.  They  arc 
resolved  automatically  by  gmtipro2  after  detection.  This  is  a  reprint  of  the  result  prepared  for  the 
Radar  conference  [8]. 
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Figure  A.  1:  A  composite  image  of  the  Ottawa  region  with  overlaid  detections/estimations  of  mov¬ 
ing  vehicles. 
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Figure  A.2:  SAR  image  of  the  Gibraltar  region  with  overlaid  ship  detections/estimations. 


DRDC  Ottawa  CR  2013-089 


25 


This  page  intentionally  left  blank. 


26 


DRDC  Ottawa  CR  2013-089 


List  of  acronyms 


ACCC  Average  Cross  Correlation  Coefficient 

ATI  Along-Track  Interferometry 

CFAR  Constant  False  Alarm  Rate 

DC  Doppler  Centroid 

DOF  Degrees  of  Freedom 

DPCA  Displaced  Phase  Centre  Antenna 

DRDC  Defence  Research  and  Development  Canada 

ECEF  Earth-Centered,  Earth-Fixed 

ENL  Equivalent  Number  of  Looks 

GLRT  Generalized  Likelihood  Ratio  Test 

LRT  Likelihood  Ratio  Test 

ML  Maximum  Likelihood 

MODEX  Moving  Object  Detection  Experiment 

MODEX1  Moving  Object  Detection  Experiment  1 

MODEX2  Moving  Object  Detection  Experiment  2 

NP  Notch  Periodogram 

PD  Probability  of  Detection 

pdf  probability  density  function 

PFA  Probability  of  False  Alarms 

PRF  Pulse  Repetition  Frequency 

PRI  Pulse  Repetition  Interval 

SAR  Synthetic  Aperture  Radar 

SMTI  Surface  Moving  Target  Indicator 

VSAR  Velocity  Synthetic  Aperture  Radar 

XML  extensible  Markup  Language 


DRDC  Ottawa  CR  2013-089 


27 


This  page  intentionally  left  blank. 


28 


DRDC  Ottawa  CR  2013-089 


3. 


DOCUMENT  CONTROL  DATA 

(Security  markings  for  the  title,  abstract  and  indexing  annotation  must  be  entered  when  the  document  is  Classified  or  Designated.) 


ORIGINATOR  (The  name  and  address  of  the  organization  preparing  the 
document.  Organizations  for  whom  the  document  was  prepared,  e.g.  Centre 
sponsoring  a  contractor’s  report,  or  tasking  agency,  are  entered  in  section  8.) 


2a. 


SECURITY  MARKING  (Overall  security 
marking  of  the  document,  including 
supplemental  markings  if  applicable.) 


UNCLASSIFIED 

2b.  CONTROLLED  GOODS 

(NON-CONTROLLED 
GOODS) 

DMC  A 

REVIEW:  GCEC  APRIL  2011 

TITLE  (The  complete  document  title  as  indicated  on  the  title  page.  Its  classification  should  be  indicated  by  the  appropriate 
abbreviation  (S,  C  or  U)  in  parentheses  after  the  title.) 


MDA  Systems  Ltd. 

75  Albert  Street,  Suite  412 
Ottawa,  ON  K1 P  5E7 


A  Class  of  CFAR  Detectors  Implemented  in  the  SAR-GMTI  Processor  gmtipro2  : 
Mathematical  Formulation  of  the  Algorithms 


4.  AUTHORS  (Last  name,  followed  by  initials  -  ranks,  titles,  etc.  not  to  be  used.) 

Dragosevic,  M.;  Burwash,  W. 


5.  DATE  OF  PUBLICATION  (Month  and  year  of  publication  of 
document.) 

6a.  NO.  OF  PAGES  (Total 
containing  information. 

Include  Annexes, 

Appendices,  etc.) 

6b.  NO.  OF  REFS  (Total 
cited  in  document.) 

February  2015 

40 

14 

7.  DESCRIPTIVE  NOTES  (The  category  of  the  document,  e.g.  technical  report,  technical  note  or  memorandum.  If  appropriate,  enter 
the  type  of  report,  e.g.  interim,  progress,  summary,  annual  or  final.  Give  the  inclusive  dates  when  a  specific  reporting  period  is 
covered.) 


Contract  Report 

8.  SPONSORING  ACTIVITY  (The  name  of  the  department  project  office  or  laboratory  sponsoring  the  research  and  development  - 
include  address.) 

Defence  Research  and  Development  Canada  -  Ottawa 
3701  Carling  Avenue,  Ottawa  ON  K1 A  0Z4,  Canada 


9a.  PROJECT  OR  GRANT  NO.  (If  appropriate,  the  applicable 
research  and  development  project  or  grant  number  under 
which  the  document  was  written.  Please  specify  whether 
project  or  grant.) 

9b.  CONTRACT  NO.  (If  appropriate,  the  applicable  number  under 
which  the  document  was  written.) 

05eq 

W771 4-0811 04/001 /SV 

10a.  ORIGINATOR’S  DOCUMENT  NUMBER  (The  official 

document  number  by  which  the  document  is  identified  by  the 
originating  activity.  This  number  must  be  unique  to  this 
document.) 

10b.  OTHER  DOCUMENT  NO(s).  (Any  other  numbers  which  may 
be  assigned  this  document  either  by  the  originator  or  by  the 
sponsor.) 

DRDC  Ottawa  CR  2013-089 

1 1 .  DOCUMENT  AVAILABILITY  (Any  limitations  on  further  dissemination  of  the  document,  other  than  those  imposed  by  security 
classification.) 


(X)  Unlimited  distribution 

(  )  Defence  departments  and  defence  contractors;  further  distribution  only  as  approved 
(  )  Defence  departments  and  Canadian  defence  contractors;  further  distribution  only  as  approved 
(  )  Government  departments  and  agencies;  further  distribution  only  as  approved 
(  )  Defence  departments;  further  distribution  only  as  approved 
(  )  Other  (please  specify): 


12.  DOCUMENT  ANNOUNCEMENT  (Any  limitation  to  the  bibliographic  announcement  of  this  document.  This  will  normally  correspond 
to  the  Document  Availability  (11).  However,  where  further  distribution  (beyond  the  audience  specified  in  (11))  is  possible,  a  wider 
announcement  audience  may  be  selected.) 

Unlimited 


13.  ABSTRACT  (A  brief  and  factual  summary  of  the  document.  It  may  also  appear  elsewhere  in  the  body  of  the  document  itself.  It  is  highly 
desirable  that  the  abstract  of  classified  documents  be  unclassified.  Each  paragraph  of  the  abstract  shall  begin  with  an  indication  of  the 
security  classification  of  the  information  in  the  paragraph  (unless  the  document  itself  is  unclassified)  represented  as  (S),  (C),  or  (U).  It  is 
not  necessary  to  include  here  abstracts  in  both  official  languages  unless  the  text  is  bilingual.) 

This  document  provides  a  detailed  mathematical  description  of  several  detection  algorithms  that 
have  been  implemented  in  the  SMTI  processor  gmtipro2  and  tested  on  the  actual  SAR  data  from 
RADARSAT-2  acquired  in  the  MODEX  mode.  The  common  feature  of  the  described  algorithms  is 
their  CFAR  property.  All  of  the  detectors  are  generalized  with  respect  to  the  number  of  available 
channels.  Thus,  they  are  defined  both  for  MODEX1  and  MODEX2  modes,  which  have  two  phys¬ 
ical  and  four  virtual  channels,  respectively.  All  of  the  detectors  work  in  the  image  domain,  but 
can  be  combined  with  DC  offset  focusing  in  order  to  integrate  all  signal  energy  and  to  improve 
the  probability  of  detection.  DC  estimation,  as  one  of  the  most  important  auxiliary  algorithms,  is 
also  presented  in  the  form  of  mathematical  equations.  The  emphasis  of  the  report  is  entirely  on 
the  mathematical  formalism  behind  the  computer  code  that  implements  these  algorithms.  Open 
literature  references  are  provided  regarding  the  derivation  and  analysis  of  these  algorithms. 


1 4.  KEYWORDS,  DESCRIPTORS  or  IDENTIFIERS  (Technically  meaningful  terms  or  short  phrases  that  characterize  a  document  and  could 
be  helpful  in  cataloguing  the  document.  They  should  be  selected  so  that  no  security  classification  is  required.  Identifiers,  such  as 
equipment  model  designation,  trade  name,  military  project  code  name,  geographic  location  may  also  be  included.  If  possible  keywords 
should  be  selected  from  a  published  thesaurus,  e.g.  Thesaurus  of  Engineering  and  Scientific  Terms  (TEST)  and  that  thesaurus  identified. 
If  it  is  not  possible  to  select  indexing  terms  which  are  Unclassified,  the  classification  of  each  should  be  indicated  as  with  the  title.) 

Ground  Moving  Target  Indication 
Synthetic  Aperture  Radar 
SAR 

C-Band  Spaceborne 

Detection 

Estimation 

Constant  False  Alarm 
CFAR 

Test  Statistic 

Generalized  Likelihood  Ratio  Test 
GLRT 

Maximu  Likelihood 
ML 

Signal  Processing 


